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SUMMARY 

A  set  of  constitutive  equations  has  recently  been  developed  to  model 
non-linear  material  behaviour  of  high  temperature  superalloys.  This  paper 
summarises  the  governing  equations  of  this  new  model  and  presents  a  simple 
FORTRAN  program  which  implements  these  equations  for  the  one¬ 
dimensional  case. 
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1. 


INTRODUCTION 


Classical  techniques  of  modelling  creep  and  plasticity  have  a  high  degree  of 
accuracy  provided  the  inelastic  strains  are  kept  small.  However,  they  are  poor 
in  describing  material  behaviour  where  significant  inelastic  strains  or  unloading 
are  involved. 

In  a  recent  investigation  [l],  a  number  of  “unified"  constitutive  models 
were  examined.  In  these  approaches,  a  single  inelastic  strain  tensor  is  used 
to  model  both  creep  and  plasticity.  A  new  constitutive  equation  was  then 
proposed. 

This  report  describes  this  new  approach,  and  discusses  a  FORTRAN  pro¬ 
gram  *  which  uses  these  new  equations  to  predict  material  behaviour.  Given 
either  stress  or  strain  values  as  input  data,  it  evaluates  the  “state  variables” , 

1. e.  drag  stress,  back  stress  etc.,  the  inelastic  strain  rate  and  the  inelastic  work. 

2.  NEW  UNIFIED  CONSTITUTIVE  MODEL 


2.1  Theory 

The  unified  constitutive  model  developed  by  Ramaswarmy  j  based  on 
the  generic  back  stress  and  drag  stress  model  proposed  by  B.  [2]  and 
Stouffer  [3].  The  Bodner-Stouffer  flow  law  can  be  written  as 


4  =  Mij  (1) 

where  A  is  a  scalar  material  function  which  can  be  found  using  a  procedure 
suggested  by  Bodner  and  Stouffer[3].  Including  a  back  stress  in  the  model 
*  Developed  by  Prof.  D.C.  Stouffer  and  V.S.  Battachar  at  the  University  of 
Cincinatti,  Ohio,  U.S.A. 
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allows  the  inelastic  strain  rate  vector  and  the  deviatoric  stress  vector  S^- 
to  have  different  directions,  and  results  in  a  flow  rule  of  the  form 


(2) 


(3) 


k2  =  \  (%  -  no)  (5o  “  n«j)  W 

as  the  second  invariants  of  the  inelastic  strain  rate  tensor  and  the  effective 
deviatoric  stress  tensor  respectively,  then  the  flow  equation  can  be  rewritten 
as 


where  D,  A  and  n  are  temperature  dependent  material  parameters.  In  this 
case  the  general  flow  equation  can  be  written  as 


This  is  the  form  of  the  flow  equation  used  in  the  FORTRAN  program. 
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The  back  stress  0  and  the  drag  stress  Z  used  in  these  equations  are  “state 
variables”  which  cannot  be  measured  directly.  They  are  represented  in  the 
form  of  evolution  equations,  where  their  rate  of  change  with  respect  to  time 
depends  upon  the  current  state  of  stress,  inelastic  strain  rate,  temperature  and 
internal  variables,  viz: 


n  =  F1(alj,i!J,nij,Z,T)  (8) 

and 

Z  =  F2(aij)efi,aiy,Z,T).  (9) 

The  multiaxial  form  of  the  back  stress  evolution  equation  used  in  the  FOR¬ 
TRAN  program  is 

where  f\,  G  and  flmal  are  scalar  parameters  determined  from  uniaxial  experi¬ 
ments,  and  where,  {or  a  uniaxial  loading  condition,  the  effective  inelastic  strain 
rate  is  defined  as 


®«//  t11) 

Also, 

fyj  =  (12) 


The  drag  stress  Z  is  a  scalar  quantity,  and  is  represented  by  the  evolution 
equation 


(13) 


Z  =  m(Z1  -  Zo)W!  +  R%. 

If  temperature  effects  are  insignificant  then  eqn  (13)  can  be  integrated  to  give 

£  =  Zi  +  (Z0  -  Zj)  e~mWl  (14) 

where  Zq  is  the  initial  value  of  Z,  Z\  is  the  saturated  value  of  drag  stress,  m 
determines  the  rate  of  cyclic  softening,  and  is  the  inelastic  work.  Further 
details  of  the  above  equations  and  their  derivations  are  given  in  reference  [l], 
whilst  reviews  of  current  methods  for  predicting  material  behaviour  are  given 
in  references  [1,4]. 

2.2  Program  Summary 

A  simple  FORTRAN  program,  which  we  will  refer  to  as  MODEL2,  has  been 
written  to  implement  this  system  of  constitutive  equations,  though  at  present 
it  is  only  applicable  to  the  one-dimensional  case.  It  may  be  run  in  either 
Stress  or  Strain  control,  with  either  isothermal  or  non-isothermal  conditions. 
The  loading  may  be  either  tensile,  torsional  or  a  combination  of  tension  or 
torsion.  Input  data  consists  of  control  variables,  load  history  points  (stress  or 
strain),  and  material  data.  A  complete  listing  of  the  input  required  can  be 
found  in  the  Appendix. 

The  subroutines  called  in  this  program,  and  a  brief  description  of  their 
functions,  are  given  in  Figure  1.  A  flow  chart  for  the  program  is  presented  in 
Figure  2. 


2.3  Detailed  Description 

•  Subroutines  INPUT1,  INPUT2  and  LOADHIST  read  in  and  initialise  ma¬ 
terial  constants,  control  parameters  and  initial  values  of  stress  or  strain. 
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•  Subroutine  TSTEP  then  determines  the  time  step  At  for  the  inelastic 
constitutive  equation,  using  current  stress  or  strain  values,  the  current  time 
step  and  control  parameters. 

The  new  time  step  At*  is  taken  as  the  minimum  of 

(i)  (At  *  A tinit)  /IE  -  10 

(ii)  (At  *  AcmaI)  /  \omax  (t  )  ~  amax  (^)l 

(iii)  (At  *  A Emai)  /  | emax  (f*)  ”  Emai  (t)j 

(iv)  A trnax  *  At, 

where  At  is  current  time  step,  Atmax  and  A tmtn  are  the  maximum  and 
minimum  time  step  factors  respectively,  At,mt  is  the  initial  time  step  and 
t*  =  t  +  At. 

However,  if  At*  <  A tm^n  *  At  then  At*  =  Atrn(n  *  At. 

•  Subroutine  INTTIME  linearly  interpolates  the  stress/strain  and  temper¬ 
atures  with  respect  to  time.  Using  subroutine  INTTEMP,  material  constants 
are  also  interpolated  with  respect  to  time  for  the  non-isothermal  case.  Sub¬ 
routine  DCCALCS  then  interpolates  the  material  stiffness  matrix  [D]  and  the 
material  compliance  matrix  (Cj. 

•  If  the  program  is  running  under  strain  control,  inelastic  strain  is  estimated 
using  the  equation 


*/;(**)  =  e/yW  + 


(15) 


•  Subroutine  STREPS,  running  under  strain  control,  calculates  elastic  strains 
and  evaluates  the  stresses  and  total  strains  as  follows 
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M0=«i/(**W/(0 


°>j(n  =  pm**)  (is) 

ei](t*)  =  ee(t*)+ej]{t‘). 

If  the  program  is  running  under  stress  control,  STREPS  calculates  only 
the  elastic  strains,  viz: 

*.(<*)  =  [<%,•(**).  (17) 

and  the  deviatoric  stresses  5,y(<*)  are  calculated  in  subroutine  DEVSTR. 

•  The  state  variables  back  stress,  drag  stress  and  K 2  are  evaluated  next, 
along  with  the  total  inelastic  strain  and  the  inelastic  work.  These  calcula¬ 
tions  are  performed  using  subroutines  DRAG,  EFFRAT,  BACK,  K2CALC, 
TOTISTR  and  WORK. 

Subroutine  DRAG  evaluates  the  drag  stress  using  eqn(14),  viz: 

Z(f)  =  Z1  +  (Z0-Z1)  e-mW‘. 

Subroutine  EFFRAT  evaluates  the  effective  inelastic  strain  rate  using 
eqn(ll),  viz: 

diff = 

Subroutine  PACK  evaluates  the  inelastic  back  stress  rate  using  eqn(10), 
and  from  this  finds  the  total  back  stress  using  eqn(l2),  viz: 
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T 


n,y(0  =  n/y(0  +  g%(0, 


(12) 


where 


«&•(**)  =  n/y(f)  +  f  (n £•(*)  +  n£(n) 


(18) 


and 


=  f\Hj  -  \i  lyr2-^/ 

^  J*max 


(10) 


Subroutine  K2CALC  evaluates  i^2i  the  second  invariant  of  the  effective 
deviatoric  stress,  from  eqn(4),  viz: 

*2  =  ^(50(t*)-ntJ(t‘))(5, :y(t*)  -niy(t*)). 

Subroutine  FLOW  evaluates  the  inelastic  strain  rate  using  eqn(7),  viz: 


«/»(**)  =  °exP 


_A  /  \n 

2  1,3*2  / 


(stJ(n  -  n,-, •(**)) 


The  total  inelastic  strain  is  found  in  subroutine  TOTISTR 


4(1‘)=4w  +  y(#)+4((‘)) 


and  the  inelastic  work  is  then  evaluated  in  subroutine  WORK 


W'  (t*)  =  W1  (t)  +  (5,y  («*)  +  St]  (())  (4(t‘)  -  e/7(t))  . 

•  Subroutine  CONCHK  checks  the  convergence  of  inelastic  strain  and  stress. 
Convergence  is  taken  to  be  achieved  if  either 


7 
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4('*) = 0 

or 

«,■(*')  =  ° 

or  if  the  calculated  percentage  error  is  less  than  the  allowable  percentage  error 
for  either  stress  or  strain. 

If  convergence  is  not  ach'jved  the  values  are  reset  in  subroutine  RESET1, 
and  the  iteration  loop  re-entered  as  shown  in  the  flow  chart,  Figure  2. 

If  convergence  is  achieved,  the  total  strain  is  then  determined  and  the 
appropriate  values  output  in  subroutine  PRINT.  Parameters  such  as  the  state 
variables  are  output  to  the  file  OUTPUT.DAT,  and  the  values  of  the  stresses 
and  strains  are  output  to  the  plot  files  PLOTl.DAT  and  PL0T2.DAT.  If  the 
analysis  is  not  yet  complete,  a  new  time  step  At  is  determined  and  the  solution 
is  recalculated. 


3.  DISCUSSION 


The  program  outlined  in  this  report  predicts  non-linear  material  behaviour 
in  the  simple  one-dimensional  case.  Predictions  using  this  model  have  been 
shown  [1,2,3]  to  be  in  good  agreement  with  experimental  results.  The  gov¬ 
erning  equations  of  the  Bodner-Stouffer  model  summarised  in  this  report  are 
applicable  to  the  three-dimensional  case.  This  program  could  be  extended  to 
handle  multiaxial  loading  conditions,  and  could  also  be  incorporated  into  a 
general  purpose  finite  element  code. 


8 


REFERENCES 


Ramaswarmy,  V.G.,  ‘A  Constitutive  Model  for  the  Inelastic  Multiaxial 
Cyclic  Response  of  a  Nickel  Base  Superalloy  Rene  80’,  NASA  Contrac¬ 
tor  Report  3998,  July  (1986). 

Bodner,  S.R.  and  Partom,  Y.,  ‘Constitutive  Equations  for  Elastic- 
Viscoplastic  Strain-Hardening  Materials’,  Journal  of  Applied  Mechan¬ 
ics,  385-389,  June  (1975). 

Stouffer,  D.C.  and  Bodner,  S.R.,  ‘A  Relationship  Between  Theory  and 
Experiment  for  a  State  Variable  Constitutive  Equation’,  American  So¬ 
ciety  for  Testing  and  Materials  STP,  765,  239-250,  1982. 

Noor,  A.K.  and  Atluri,  S.N.,  ‘Advances  and  Trends  in  Computational 
Structural  Mechanics’,  AIAA  Journal,  25,  7,  977-995,  1987. 


SUBROUTINES : 


INPUT1 

INPUT2 

LOADHIST 


I NTT I ME 
I NTTEMP 
DCCALC 
TSTEP 
STREPS 


DEVS  TR 
DRAG 

EFFRAT 

BACK 

K2CALC 

FLOW 

TOTISTR 

WORK 

CONCHK 

RESET1 

RESET2 

ERRPRT 

PRINT 


F igure  1  : 


-  Reads  Control  Data 

-  Reads  Material  Constants 

-  From  control  data  this  sets  up  the 
stress/strain  and  temperature  variation 
curve  wrt  time 

-  Interpolates  stress/strain  &  temperature 
wrt  time 

-  Interpolates  material  constants  wrt 

temperature 

-  Calculates  D  and  C  matrix 

■  Calculates  current  -  new  time  step  dt 

-  Calculates  stresses  if  strain  control 
and  elastic  strains  if  stress  control 

-  Calculates  deviatoric  stresses  at  tdi 

-  Calculates  Drag  stress  7.  at.  tdt 

-  Calculates  effective  inelastic  strain 
rate  at  tdt 

-  Calculates  Back  stress  Omega  at  tdt 

-  Calculates  value  of  AK2  at  tdt 

-  Calculates  inelastic  strain  rate  at  tdt 

-  Calculates  total  inelastic  strain  at  tdt 

-  Calculates  inelastic  work  at  tdt. 

-  Check  Convergence  of  values 

-  Resets  values  if  there  is  no  convorgei.ce 

-  Resets  values  if  there  is  convergence 

-  Prints  error  messages  into  a  file 

-  Prints  output  into  a  file 


Subroutines  called  in  MODEL2 


INPUT  DATA. 
DETERMINE  LOAD 
HISTORY 


NOTE:  Time  t*=t+At 


STOP 


FIGURE  2.  FLOW  CHART  FOR  PROGRAM  MODEL  2.F 


APPENDIX:  Variables  used  in  the  program  M0DEL2 : 


INPUT  VARIABLES: 

1.  TITLES: 

*ETITLE  -  title  line  in  the  file  INPUT . DAT 
“MTITLE  -  title  line  in  the  file  MATL.DAT 

2.  CONTROL  VARIABLES: - 

*NCTL  -  decides  stress/strain  control  . 

=  1  stress  control 
=  2  strain  control 
‘NTYPE-  decides  the  load  type 
=  1  tensile  load  only 
=  2  torsional  load  only 
=  3  tension  +  torsion 

‘NTHERM-  decides  isothermal/r.onisotherma  1  case 
=  1  isothermal 
=  2  nonisothermal 

‘NCYTEN-  number  of  cycles  in  tension  loading  =  or  >1 
*NCYTOR-  number  of  cycles  in  torsion  loading  =  or  >1 
*NCYTEM-  number  of  cycles  in  temperature  loading 
=  1  isothermal 

>  or  =  1  nonisothermal 

*NPTEN  -  number  of  points  in  tension  loading  in 
first  cycle  =  or  >1 

*NPTOR  -  number  of  cycles  in  torsion  loading  in 
first  cycle  =  or  >1 

•NPTEM  -  number  of  cycles  in  temperature  loading  in 
first  cycle  =  or  >1 

LASTEN\ 

LASTOR  >  last  points  in  respective  loadings 
LASTEM/ 

*TSTOP  -  time  at  which  analysis  is  stopped  =  or  >0 
‘EST0P1-  strain  in  the  tension  load  direction  at  which 
analysis  is  to  be  stopped  =  or  >0 
‘EST0P2-  strain  in  the  shear  stress  direction  at  which 
analysis  is  to  be  stopped  =  or  >0 

*DTPR  -  output  printed  into  OUTPUT.DAT  once  in  DTPR  time 
steps  (time  conli oiled  output) 

=  0  if  not  time  controlled 

>  O  if  time  controlled 

*DEPR1  -  output  printed  once  in  DEPR1  increase  in  strain 

in  the  tension  direction  (strain  controlled  output) 
=  0  if  not  stra in - i n- ax i a  1  -  d  irect  ion  controlled 

>  0  if  stra in  -  in  -  ax ia 1 -d irect ion  controlled 
-  output  printed  once  in  DEPR2  increase  in  strain 

in  the  shear  direction  (strain  controlled  output) 

=  O  if  not  strain  in  shear -d i rect i on  controlled 

>  O  if  stra in  -  in  -  shear - d irect i on  controlled 


*DEPR2 


*CTR (20) -control  variables 

(1)  -maximum  stress  step  allowable  in  analysis 

(2)  -maximum  inelastic  strain  step 

(3)  -maximum  inelastic  strain  rate  step 

(4)  -maximum  time  step  factor 

(5)  -minimum  time  step  factor 

(6)  -initial  time  step 

CTR(l-6)  values  >  0 


3.  ERROR  CONTROL  DATA: 

*PERSTR  -allowable  %  error  in  stress  calculations 
*PEREPS  -allowable  %  error  in  strain  calculations 
‘MAXITER-maximum  number  of  iterations  allowable  for 
convergence 

4.  MATERIAL  CONSTANTS  DATA: 

“NUMCON  -number  of  input  material  constants (=  or  <20) 
**NUMTEM  -number  of  temperatures  at  which  material 
constants  are  being  input (=  or  <10) 

* ‘CONS (20, 10) -material  constants  . 

CONS ( I , J)  I  =  1,20  J  =  1,10 

1  =  1  :  A 

2  :  D 

3  :  YOUNGS  MODULUS 

4  :  fl 

5  :  g 

6  :  m 

7  :  n 

8  :  poissons  ratio 

9  :  ZO 

10  :  Z1 

11  :  OMEGA  MAX . 

CONINT ( 20) - interpo lated  material  constants 

D(6,6)  -material  stiffness  matrix 
C(6,6)  -material  compliance  matrix 


5.  INPUT  LOAD  DATA: (strain  inputs  are  tensor  strain  values) 

‘DXTTEN ( 10 , 2)  -  tensile  load  points  in  the  first  cycle 

(stress/strain) 

‘DXTTOR ( 10, 2 )  -  torsion  load  points  in  the  first  cycle 

(stress/strain) 

‘DXTTEM ( 10, 2 )  -  temperature  load  points  in  the  first  cycle 

DXT . . . (1,1)  PARAMETER  VALUE 
DXT . ..(1,2)  TIME  VALUE 


6.  TEMPERATURE  VARIABLES: 


** TEMDAT ( 10) -  temperatures  at  which  material  constants  have 
been  input 

NOTE:  1.  ALL  STRAIN  INPUTS  are  TENSOR  values  except  the  values 

of  the  strains  which  are  printed  into  PL0T1  and 
PL0T2.DAT  files;  these  are  engineering  strain  values. 

2.  A  before  the  variable  indicates  it  is  an  input 

from  the  file  'INPUT.DAT' . 

3.  A  '**'  before  the  variable  indicates  it  is  an  input 

from  the  file  'MATL.DAT*. 


INTERNAL  AND  OUTPUT  VARIABLES: 


7.  CALCULATED  LOAD  DATA: 

DATTEN (500, 2)  -  data  points  in  the  entire  load  history 
DATTOR (500, 2)  -  in  the  respective  cycles. 

DATI EM ( 500 , 2 )  - 

DAT ...(1,1)  PARAMETER  VALUE 
DAT.  ..(1.2)  TIME  VALUE 

8.  STRESS  VARIABLES: 

STR_t(6)  -stress  at  time  t 
STR_tdt (6)  -  stress  at  time  t  +  dt 

PSTR(6)  -stress  at  t+dt  stored  from  the  last  iteration 
DEVS_t (6)  -deviatoric  stress  at  time  t 
DEVS_tdt (6)  -deviatoric  stress  at  time  t+dt 

9.  STRAIN  VARIABLES: (all  strains  are  tensor  strains  !) 

EPS_t (6)  -total  strain  at  time  t 
EPS_tdt (6) - tota 1  strain  at  time  t+dt 

EPSE_t (6) -elastic  strain  at  time  t 
EPS_tdt (6) -elastic  strain  at  time  t+dt 

EPSI_t(6)  -inelastic  strain  at  time  t 
EPSI_tdt (6) -  inelastic  strain  at  time  t  +  dt 
PEPSI_tdt (6) -  inelastic  strain  at  time  t  +  dt  stored  from 
the  last  iteration 

EPSID_t (6) - inel ast ic  strain  rate  at  time  t 
EPSID_tdt (6)  -  inelastic  strain  rate  at  time  t  +  dt 
PEPSID (6) -  inelastic  strain  rate  at  time  t  +  dt  stored 
from  the  last  iteration 

EEIEPSD-e f feet ive  inelastic  strain  rate  at  time  t+dt 


10. STATE  VARIABLES: 


Z_t  -drag  stress  at  time  t 
Z_tdt  -drag  stress  at  time  t+dt 

0M__t  (6)  -total  back  stress  at  time  t 
0M_tdt (6)  -total  back  stress  at  time  t+dt 

POM (6)  -total  back  stress  at  time  t+dt  from  the  previous 
iteration 

0MI_t(6)  -inelastic  back  stress  at  time  t 
0MI_tdt (6) - inelastic  back  stress  at  time  t+dt 
POMI (6)  -inelastic  beck  stress  rate  at  time  t+dt  from  the 
previous  iteration 

0MID_t (6) - inelastic  back  stress  rate  at  time  t 
0MID_tdt (6) - inelastic  back  stress  rate  at  time  t+dt 
P0MID(6)  -inelastic  back  stress  rate  at  time  t+dt  from 
the  previous  iteration 

AK2  -  K2  value  at  time  t+dt 

11.  TIME  VARIABLES: 

TIME  -time  for  the  present  time  step  (same  as  t+dt) 
dt  -time  increment  (=TIME-previous  time) 

TPRINT-ti me  at  which  output  was  printed  last 

12.  WORK  VARIABLES: 

WRK_t  -total  inelastic  work,  at  time  t 

WRK_tdt  -total  inelastic  work  at  time  t+dt 

DELWRK  -increment  of  work  from  time  t  to  time  t+dt 

13.  TEMPERATURE  VARIABLES: 

TEMP  * temperature  at  time  t+dt 

TEM1,TEM2  -two  input  temperatures  between  which  TEMP 

lies  (used  for  interpolation  of  matl  const.) 


14. OTHER  VARIABLES: 

ITER  -number  of  iterations 
NCODE-convergence  code 
=0  converged 
=1  not  converged 

IERCOD-code  passed  on  to  error  printing  subroutine 
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